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In this talk we discuss an algorithm for the numerical calculation of one-loop QCD amplitudes 
and present results at next-to-leading order for jet observables in electron-positron annihilation 
calculated with the above-mentioned method. The algorithm consists of subtraction terms, ap- 
proximating the soft, collinear and ultraviolet divergences of QCD one-loop amplitudes, as well 
as a method to deform the integration contour for the loop integration into the complex plane 
to match Feynman's i8 rule. The algorithm is formulated at the amplitude level and does not 
rely on Feynman graphs. Therefore all ingredients of the algorithm can be calculated efficiently 
using recurrence relations. The application of this method to the leading-colour contribution of 
e + e~ — > n jets, with n up to seven, demonstrates the efficiency of the approach. 
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1. Introduction 



The development of a fully numerical algorithm to calculate multi-parton QCD amplitudes 
and collider observables at next-to-leading order (NLO) accuracy and the implementation of such 
an algorithm in a viable Monte Carlo matrix element generator is a rather involved subject. How- 
ever, the need for generators of this kind is motivated by the need for accurate QCD background 
predictions to LHC physics. On inclusion of further Standard Model interactions such a generator 
can also be used to directly compute other relevant processes in addition to QCD background. Our 
focus lies upon the virtual part of the NLO calculation, i.e. on evaluating the one-loop integration 
numerically, where we employ and extend the ideas of the subtraction method to the virtual part 
[]|-|tl- I n this regard our algorithm is different from the commonly used approaches, based on cut 
techniques and generalised unitarity or on more traditional Feynman graph approaches [^-[T^], but 
shows promising features especially for the implementation in a numerical program. The algorithm 
consists of local subtraction terms to subtract divergences arising from the soft, collinear and ultra- 
violet regions of the virtual part, which render the integrand finite in the respective regions, and of a 
method to deform the integration contour of the loop integration into the complex plane in order to 
circumvent the remaining on-shell singularities. It works on the level of colour-ordered primitive 
amplitudes, where we utilise recursive algorithms to compute the corresponding one-loop off-shell 
currents for the bare primitive amplitudes, and is therefore fast and easily implemented. 

The local subtraction terms for the soft and the collinear regions are formulated directly on the 
amplitude level. These subtraction terms are proportional to the corresponding Born amplitudes 
and are easily implemented. The local subtraction terms for the ultraviolet region are known to in- 
clude only propagator and vertex corrections, where the corresponding graphs are expanded around 
a new ultraviolet propagator. The total local ultraviolet subtraction term is constructed recursively 
from the propagator and vertex subtraction terms. Once the local subtraction terms are applied to 
the bare integrand it can be integrated numerically in four dimensions in loop-momentum space. 
However, singularities still remain, since one or more of the propagators still may go on-shell for 
certain real values of the loop-momentum. To avoid these singularities we deform the integration 
contour into the complex space. The contour deformation can be implemented in two ways, either 
by a direct deformation of the loop four-momentum only or alternatively by introducing Feyn- 
man parameters and deforming the loop four-momentum as well as the corresponding Feynman 
parameters. The numerical loop integration is performed together with the integration over the 
phase-space of the external particles in one Monte Carlo integration. The subtraction terms yield 
simple results upon analytic integration over the loop-momentum, and the resulting pole structures 
cancel exactly against the pole structures from the soft and collinear parts of the real emission 
contributions as well as of the ultraviolet countertem from renormalisation. Hence, the algorithm 
goes hand in hand with the usual subtraction method for the real emission contributions, where we 
employ Catani-Seymour dipole subtraction [20-24]. As a proof of concept the algorithm has been 
tested so far for e + e~ — > n jets, with n up to seven, for the massless case in the large-A^. limit. Up 
to four jets we reproduce the known results for the respective jet rates with very good agreement. 
Increasing the number of jets up to seven shows a good scaling behaviour in CPU time with respect 
to the number of final state particles. 
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2. Subtraction method 

The subtraction method is widely used to render the real emission part of a NLO calculation 
suitable for a numerical Monte Carlo integration. The contributions to an infrared-safe observable 
at next-to-leading order with n final state particles can be written as 

(0) NL0 = J O n+l do R + J O n do v + J O n do c . (2.1) 

n+l n n 

Here a rather condensed notation is used. do R denotes the real emission contribution, whose matrix 
elements are given by the square of the Born amplitudes with [n + 3) partons |A^ 3 1 2 . da v denotes 
the virtual contribution, whose matrix elements are given by the interference term of the one-loop 
and Born amplitude 3i(A^ 2 A^ 2 ) and da c denotes a collinear subtraction term, which subtracts 
the initial state collinear singularities. Each term is separately divergent and only their sum is finite. 
One adds and subtracts a suitably chosen piece to be able to perform the phase space integrations 
by Monte Carlo methods: 



(0) NL0 = J {O n+l da R - O n do A ) + J O n do v + O n do c + O n J do A \ . < 2.2 ) 

n+l n \ 1 



The first term [O n+ \do R — O n da A ) is by construction integrable over the (n + l)-particle phase 
space and can be evaluated numerically. After integration of the subtraction term over the unre- 
solved one-parton phase space the infrared divergences of the virtual contribution from the one-loop 
amplitude cancel with the infrared poles of the subtraction terms. Therefore the second term is also 
infrared finite and can be evaluated numerically, provided the analytical result of the one-loop 
amplitude is known. 

We extend this subtraction method to the virtual part such that we can evaluate the one-loop 
integral of the one-loop amplitude numerically. The renormalised one-loop amplitude is related to 
the bare amplitude by 

^ (1) = <e + ^cx\ (2-3) 

where denotes the ultraviolet counterterm from renormalisation. The bare amplitude involves 
the loop integration 

(1 ) _ [_dPk_yM 24 

bare ~ J (2n) D bare ' ( j 



where denotes the integrand of the bare one-loop amplitude. We introduce subtraction terms 



which match locally the singular behaviour of the bare integrand: 

^(D _ f (yw 

-^bare + -^CT — J (2n) D V bare soft co11 UV / 
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Analogous to ^are' tne integrands of the subtraction terms si^ are denoted by \ where x 
is equal to soft, coll or UV. The expression in the first bracket is finite and can therefore be 
integrated numerically in four dimensions. The integrated subtraction terms in the second bracket 
can be easily calculated analytically in D dimensions. The poles in the dimensional regularisation 
parameter of the integrated subtraction terms are cancelled by the corresponding poles from the 
ultraviolet counterterms, initial state collinear subtraction terms and the integrated real emission 
subtraction terms. 

In analogy to the one-loop amplitude we can write do v = docj + J jj^do^ and then the 
NLO contributions reads 

{0) NL0 = J (O n+l do R - O n da A ) + J [o„fioLc - Onde A ') 

n+\ rc+loop 

+ j O n do C j + O n do c + O n J do A + O n J do A 

n \ loop 



In a condensed notation this reads 

(0) NW = (Off + <« a i + (Con. (2-7) 
Every single term is finite and can be evaluated numerically. 



3. Local infrared subtraction terms 

Amplitudes in QCD may be decomposed into group-theoretical factors (carrying the colour 
structures) multiplied by kinematic factors called partial amplitudes. At the loop level partial am- 
plitudes may further be decomposed into primitive amplitudes. A few important properties of 
primitive amplitudes shall be given: Firstly, primitive amplitudes are gauge invariant. This is im- 
portant for the formal proof of the method. Secondly, for a given number of external legs primitive 
amplitudes have a fixed cyclic ordering of the external legs and a definite routing of the external 
fermion lines through the loop. This ensures that each propagator in the loop is uniquely defined in 
type, be it a quark or a gluon/ghost propagator, and position. Due to the fixed cycling ordering there 
are only n different loop propagators occuring in a primitive amplitude with n external legs. With 




Pi Pn 

Figure 1: The labelling of the momenta for a primitive one-loop amplitude. 
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n 3 
the notation as in fig. |1] we define kj = k — qj, with qj = £ where k is the integration variable 

i=i 

of the one-loop integral and the p,'s are the external momenta. We can write the bare primitive 
one-loop amplitude as 

= / j2nr G ^ G ^ = m) flq-mj+is- {3A) 

P a {k) is a polynomial of degree a in the loop momenta k and the +/5-prescription tells us in which 
direction the poles of the propagators should be avoided. 

Soft singularities arise for k ~ qi, whenever pf = m,- = 0, pf +1 = mj +l . In this case we 
have a massless particle exchanged between two on-shell particles and the momentum kj is soft. 
Collinear singularities arise for k ~ qt —xpu whenever p\ = 0, m,_i = 0, m\ = 0, where x G [0, 1]. 
In this case a massless external on-shell particle is attached to two massless propagators and the 
momenta i , kj and p\ are collinear. The soft and collinear subtraction terms for massless QCD 
read 

G (i) _ 4i y Pj-Pj+i AO) 

U soft ~ 4 ' L u2 k 2 k 2 A i ' 



S jguv {kj_ v kj) Sj+i guv (kj,kj +l ) 



Af , (3.2) 



where the sum over j € I g is over all gluon propagators j inside the loop. Furthermore, Sj = 
1 if the external line j corresponds to a quark and Sj = 1/2 if it corresponds to a gluon. The 
function guv ensures that the integration over the loop momentum is ultraviolet finite. The soft and 
collinear subtraction terms are formulated directly on the amplitude level and are proportional to 
the corresponding Born amplitudes. Upon integration they yield simple analytic results: 

with S e = (47r) e exp(— BJe) the typical volume factor in dimensional regularization, where y E de- 
notes the Euler-Mascheroni constant, /I denotes the renormalization scale in dimensional regular- 
ization and £ is defined through D = 4 — 2e. 

The ultraviolet subtraction terms correspond to propagator and vertex corrections. The sub- 
traction terms are obtained by expanding the relevant loop propagators around a new ultraviolet 
propagator (k 2 — jli^ v ) , where k = k — Q: For a single propagator we have 

_JL^ = ' + g • (g - o: [Ai ~ e> 2 +if + P* : (g - a] 2 + g ( ' ) , 3 . 4) 
(t-p) 2 ' 2 -m5v (P-mSv) (* 2 -fuv) (* 2 -mSv) VWV 

We can always add finite terms to the subtraction terms. For the ultraviolet subtraction terms we 
choose the finite terms such that the finite parts of the integrated ultraviolet subtraction terms are 
independent of Q and proportional to the pole part, with the same constant of proportionality for 
all ultraviolet subtraction terms. This ensures that the sum of all integrated UV subtraction terms 
is again proportional to a tree-level amplitude. 
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4. Contour deformation 

Having a complete list of ultraviolet and infrared subtraction terms at hand, we can ensure 
that the integration over the loop momentum gives a finite result and can therefore be performed 
in four dimensions. However, this does not yet imply that we can safely integrate each of the 
four components of the loop momentum k^ 1 from minus infinity to plus infinity along the real axis. 
There is still the possibility that some of the loop propagators go on-shell for real values of the loop 
momentum. If the contour is not pinched this is harmless, as we may escape into the complex plane 
in a direction indicated by Feynman's +/5-prescription. However, it implies that the integration 
should be done over a region of real dimension 4 in the complex space C 4 . Let us consider an 
integral corresponding to a primitive one-loop amplitude with n propagators minus the appropriate 
IR- and UV-subtraction terms: 

dA ~ k 6,(1) r dA ~ k — " 1 



(2%fV^ " soft ~«* "™J J (liu^^Jcj-mj + iS (41) 

where P(k) is a polynomial of the loop momentum k 11 and the integration is over a complex contour 
in order to avoid whenever possible the poles of the propagators. We discuss the method of the 
direct deformation of the loop momentum. We set 

k = k + iK(k) (4.2) 

where W 1 is real. After this deformation our integral equals 



d 4 k 

(my 



dk» 



dk v 



P(k(k))Y\-o 5 — K • (4-3) 



To match Feynman's +/c> -prescription we have to construct the deformation vector K such that 

kj - m) = -> k r K > 0. (4.4) 

We remark that the numerical stability of the Monte Carlo integration depends strongly on the 
definition of the deformation vector K\ 



5. Recursion relations 

We use Berends-Giele type recursion relations [ ^ ] to compute the tree amplitude, the bare 
one-loop integrand G^ e an d the total UV subtraction term G^y. These recursion relations are 
shown in fig. ||] for the case of a three- valent toy model. 



6. NLO results for n-jets in electron-positron annihilation 

We have calculated results for jet observables in electron-positron annihilation, where the jets 
are defined by the Durham jet algorithm. The cross section for n jets normalised to the LO cross 
section for e + e~ — > hadrons reads 
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Figure 2: Recursive relations for a three-valent toy model. 



One can expand the perturbative coefficient A„ and B n in \/N c : 

n-2 



A n = N c 



A n ,\c + • 



1 



B„ 



N c 



n-l 



B n ,\c + 1 



(6.2) 



We calculate the leading order coefficient A n ^ c and the next-to-leading order coefficient B, h \ c for 
n < 7 at the renormalisation scale /I equal to the centre of mass energy. We take the centre of 
mass energy to be equal to the mass of the Z-boson. The scale variation can be restored from the 
renormalisation group equation. The calculation is done with five massless quark flavours. Fig. || 
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Figure 3: Comparison of the NLO corrections to the two-, three- and four-jet rate between the numerical 
calculation and an analytic calculation. The error bars from the Monte Carlo integration are shown and are 
almost invisible. 

shows the comparison of our numerical approach with the well-known results for two, three and 
four jets [^, ^7|]. We observe an excellent agreement. The results for five, six and seven jets for 
the jet parameter y cut = 0.0006 are 

^A 5 , lc = (2.4764 ±0.0002) • 10 4 , tI B 5,1c = (1.84±0.15) • 10 6 , 

8 16 

^A 6; i c = (2.874 ±0.002) • 10 5 , ^fi 6 , lc = (3.88 ±0.18) • 10 7 , 

16 51 

^A 7 ,ic = (2.49 ± 0.08) -10 6 , ^ffi 7 ,ic = (5.4±0.3)-10 8 . (6.3) 

32 o4 
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